Most modern dynamic models are constructed based on finite spatial discretizations of continuous systems, often resulting in a considerable number of degrees of freedom in the model. Consequently, for fast and efficient estimation of dynamic behavior, as well as optimizations and closed-loop control designs, a model reduction is typically performed. Desirable characteristics of a reduced-order dynamic model include that the size of the system must not be too large, the model must be robust and have a good fidelity, it must be in the state-space, time domain formulation for implementation of active control systems and nonlinear time analysis, and finally, the reduction process itself must not be too expensive.
There have been many model reduction methods available, but most of them require modifying the main frame of the computational model, and are prone to a long model construction time if the model is subjected to many driving inputs. The latter in particular is true in unsteady Computational Fluid Dynamic (CFD) applications where the moving solid boundary is often described by many structural mode inputs. For example, a typical commercial airplane simulation may involve as many as 200 structural modes.
Recently, Eigensystem Realization Algorithm (ERA) was used successfully in aeroelastic flutter predictions using a computational fluid dynamics code (CFD) developed by the NASA Langley Research Center known as CFL3D. (See Juang, J. N., Applied System Identification, Prentice Hall Englewood Cliffs, N.J., 1994; Silva, W. A., and Bartels, R. E., Development of Reduced-Order Models for Aeroelastic Analysis and Flutter Prediction Using CFL3Dv6.0 Code, AIAA-2002-1596; and Hong, M. S., Bhatia, K. G., SenGupta, G., Kim, T., Kuruvila, G., Silva, W. A., Bartels, and R., Biedron, R., Simulations of a Twin-Engine Transport Flutter Model In the Transonic Dynamics Tunnel, IFASD Paper 2003-US-44). The ERA method, which is usually used as a system identification technique, has a very attractive feature in that unlike model reduction methods based on Galerkin scheme (e.g., Dowell, E. H., Hall, K. C., and Romanowski, M. C., Eigenmode Analysis in Unsteady Aerodynamics: Reduced-Order Models, Applied Mechanics Review, Vol. 50, No. 6, 1997, pp. 371-386), there is no need for on-line implementation of the algorithm. That is, the ERA method is a post-processing tool that identifies and generates system matrices based on the input and output data alone.
FIG. 1 is a schematic view of the ERA (a.k.a. Pulse/ERA) method 100 of modeling a dynamic system in accordance with the prior art. For simplicity, only its fundamental state-space realization theory which is attributed to Ho and Kalman is discussed. (see Ho, B. L. and Kalman, R. E., Effective Construction of Linear State-Variable Models from Input/Output Data, Proceedings of the 3rd Annual Allerton Conference on Circuit and System Theory, 1965, pp. 449-459). For a general description of the Pulse/ERA method, see Juang, J.-N. and Pappa, R. S., An Eigensystem Realization Algorithm for Modal Parameter Identification and Model Reduction, Journal of Guidance, Control, and Dynamics, Vol. 8, 1985, pp. 620-627.
In the Pulse/ERA method 100, it is assumed that the dynamic system under consideration can be described in the following finite-dimensional, discrete-time form:xn+1=Axn+Bun  (1)yn=Cxn+Dun  (2)
wherex≡(L×1) state vector  (3)u≡(Ni×1) input vector  (4)y≡(N0×1) output vector  (5)The system matrices (A, B), (A, C) are assumed controllable and observable. First, given M+1 equally distributed time steps,
                              t          n                ≡                  n          ⁢                                          ⁢          Δ          ⁢                                          ⁢          t                                              (                                    n              =              0                        ,            1            ,            2            ,            …            ⁢                                                  ,            M                    )                ,            for a single i-th input vector the system output is sampled subjected to the unit pulse,
                              u          i          n                =                              δ            ni                    ≡                      {                                                            1                                                                      (                                          n                      =                      0                                        )                                                                                                0                                                                      (                                                                  ni                        =                        1                                            ,                      2                      ,                      …                      ⁢                                                                                          ,                      M                                        )                                                                        }                                              (        6        )            Assuming zero initial condition, x0=x(0)=0, one obtains
                                                                                             y                  1                  0                                =                                ⁢                                  d                  i                                                                                                                          y                  i                  1                                =                                ⁢                                  Cb                  i                                                                                                                          y                  i                  2                                =                                ⁢                                  CAb                  i                                                                                                                          y                  i                  3                                =                                ⁢                                                      CA                    2                                    ⁢                                      b                    i                                                                                                                          ⋮                ⁢                                ⁢                ⋮                                                                                                          y                  i                  M                                =                                ⁢                                                      CA                                          M                      -                      1                                                        ⁢                                      b                    i                                                                                                                                          ⁢                where                                                                                                          b                  i                                ≡                                ⁢                                  i                  ⁢                                      -                                    ⁢                  th                  ⁢                                                                          ⁢                  column                  ⁢                                                                          ⁢                  of                  ⁢                                                                          ⁢                  B                                                                                                                          d                  i                                ≡                                ⁢                                  i                  ⁢                                      -                                    ⁢                  th                  ⁢                                                                          ⁢                  column                  ⁢                                                                          ⁢                  of                  ⁢                                                                          ⁢                  D                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        (                                  7                                  )                                                                                                                                                                                                                                                                                                                                                                                                                                                                                (                      8                      )                                                                                                                                              (                9                )                                                        The constant matrices in the above sequence are known as the Markov parameters. As shown in FIG. 1, for all inputs 102 into the dynamic system, the Markov parameters are computed (block 104) creating the sequence of pulse-response matrices:
                                                                        Y                n                            ≡                                                                                                                                  b                        1                        n                                                                                                            y                        2                        n                                                                                    …                                                                                                                y                          Ni                          ni                                                ⁢                                                                                                                                                                      (                                                n                  =                  0                                ,                1                ,                2                ,                                  …                  ⁢                                                                          ⁢                  M                                            )                                                          (        10        )            Next, based on the system Markov parameters (block 104), the Pulse/ERA method 100 defines N0×(Ni×(M−1)) Hankel matrices:
                              H          0                ≡                                                                                        Y                  1                                                                              Y                  2                                                            …                                                              Y                                      M                    -                    1                                                                                                ⁢                                  ⁢                                  ≡                  C          ⁢                                                                                  I                                                  A                                                                      A                    2                                                                    …                                                                                                                    A                                                  M                          -                          2                                                                                                            ⁢                    B                                                                                                          (        11        )                                          H          1                ≡                                                      Y              2                        ⁢                                                  ⁢                          Y              3                        ⁢                                                  ⁢            …            ⁢                                                  ⁢                          Y              M                                                ⁢                                  ⁢                                  ≡                  C          ⁢                                                                                  A                                                                      A                    2                                                                    …                                                                                                                    A                                                  M                          -                          1                                                                                                            ⁢                    B                                                                                                          (        12        )            
Singular-value-decomposition (SVD) of H0 yields
                                          H            0                    ≡                      U            ⁢                                                  ⁢                          ∑                                                          ⁢                              V                T                                              ⁢                                          ≃                                                                                                                                                        U                        R                                                                                                            U                        D                                                                                                                        ⁡                              [                                                                                                    ∑                        R                                                                                    0                                                                                                  0                                                              0                                                                      ]                                      ⁡                          [                                                                                          V                      R                      T                                                                                                                                  V                      D                      T                                                                                  ]                                      ⁢                                  =                              U            R                    ⁢                                                    ∑                                  1                  /                  2                                            R                        ⁢                                                            ∑                                      1                    /                    2                                                  R                            ⁢                              V                R                T                                                                        (        13        )            where R≡rank(H0). Finally, a balanced realization of the system under question is obtained by pseudo-inverting various submatrix components (block 108 of FIG. 1) as follows:
                    D        =                              Y            0                    ⁡                      (                                          N                0                            ×              1                        )                                              (        14        )                                C        ≃                              U            R                    ⁢                                    ∑              R                              1                /                2                                      ⁢                          (                                                N                  0                                ×                R                            )                                                          (        15        )                                B        ≃                  the          ⁢                                          ⁢          first          ⁢                                          ⁢                      N            i                    ⁢                                          ⁢          columns          ⁢                                          ⁢          of          ⁢                                          ⁢                                                    ∑                                  1                  /                  2                                            R                        ⁢                                          V                R                T                            ⁢                                                          ⁢                              (                                  R                  ×                                      N                    i                                                  )                                                                        (        16        )                                A        ≃                                            ∑                                                -                  1                                /                2                                      R                    ⁢                                    U              R              T                        ⁢                          H              1                        ⁢                          V              R                        ⁢                                                            ∑                                                            -                      1                                        /                    2                                                  R                            ⁢                                                          ⁢                              (                                  R                  ×                  R                                )                                                                        (        17        )            
Since R<<L, the above model represents a reduced-order realization of the original system. Note that the realization is optimal in that it is balanced between inputs and outputs. However, the total number of samples taken is Ni×(M+1), which increases proportional to the number of inputs. Also, for an accurate system identification H0 must have sufficient columns and rows, i.e., Ni×(M−1), ≧R and N0≧R.
For a very large data set Yn with many time steps and a large number of inputs, the Eigensystem Realization Algorithm/Data Correlations (ERA/DC) method may be preferred. As described, for example, in Juang, J. N., Applied System Identification, Prentice Hall Englewood Cliffs, N.J., 1994, the ERA/DC method may be used to compress the amount of data and reduce the computation time required for the SVD of the Hankel matrix.
Although desirable results have been achieved using the prior art computational methods, there is room for improvement. For example, if the unsteady CFD model is driven by multiple structural inputs, as described above, the computation time required to obtain all the pulse responses increases proportional to the number of the inputs, making the ERA method undesirably slow and inefficient. Therefore, a need exists for improved methods for model reduction and system identification of large-scaled linear dynamic systems having multiple inputs